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To understand the role of local sublattice imbalance in low energy spectra of s = | quantum antiferromag- 
nets (as observed by [L. Wang and A.W. Sandvik, Phys. Rev. Lett. 97, 1 17204 (2006); Phys. Rev. B 81, 054417 
(2010)]), we study the s = k quantum nearest neighbor Heisenberg antiferromagnet on the coordination 3 
Cayley tree. We perform many-body calculations using an implementation of the Density Matrix Renormaliza- 
tion Group (DMRG) technique for generic tree graphs. We discover that the bond-centered Cayley tree has a 
quasi-degenerate set of low lying tower of states and an "anomalous" singlet-triplet finite size gap scaling. For 
understanding the construction of the first excited state from the many-body ground state, we consider a wave- 
function ansatz given by the single mode approximation (SMA), which yields a high overlap with the DMRG 
wavefunction. Observing the ground state entanglement spectrum helps us to analytically understand the scaling 
of the finite size spin gap and leads us to a picture of the low energy degrees of freedom being "giant spins" aris- 
ing out of sublattice imbalance. The Schwinger Boson mean field theory has been generalized to non uniform 
lattices and ground states have been found which are spatially inhomogeneous in the mean field parameters. 



I. INTRODUCTION 

Quantum antiferromagnetism for unfrustrated systems has 
been one of the most extensively researched subjects in con- 
densed matter physics. One of the simplest models in this 
family, the nearest neighbor spin i Heisenberg antiferromag- 
net on the square lattice, has been studied extensively: analyt- 
ically with Spin Wave 1-3 and Schwinger Boson approaches 4-8 
and numerically with Quantum Monte Carlo 9 (which has no 
"sign problem" for bipartite lattices). That said, effects from 
physical imperfections such as the presence of open edges 10 ' 1 1 
and static non magnetic impurities 12 are less well understood 
and hence are areas of active research. Random dilution, 
for example, has been found to affect the low energy spec- 
trum drastically, as has been demonstrated in recent numer- 
ical studies by Wang and Sandvik 13 . They showed that the 
singlet-triplet energy gap on finite clusters, at the percola- 
tion threshold, is much smaller than that expected from the 
conventional Anderson tower of states (or "quantum rotor" 
states) 14-18 . This was in disagreement with a previous theo- 
retical result 19 . 

Wang and Sandvik attributed the mechanism of the creation 
of this low energy scale to the presence of regions of "local 
sublattice imbalance", arising out of areas in the cluster where 
the number of opposite sublattice sites were unequal. In their 
picture, geometric constraints of a diluted square lattice for- 
bid spins from pairing up with their neighbors into dimers, 
leaving some of them unpaired or "dangling" 13 . They believe 
these emergent spins to be the effective low energy degrees of 
freedom. 

By considering the undiluted nearest neighbor Heisenberg 
model on the Cayley tree (or Bethe Lattice) (shown in Fig. 1) 
we have been able to systematically investigate the effect of 
sublattice imbalance on the low energy spectrum. In partic- 
ular, we contend that the effect of sublattice imbalance is to 
create a "tower of states", lower than the Anderson tower of 
states. Aided by numerical calculations, we propose a frame- 
work for understanding this relationship. We also find that 
Schwinger Boson Mean Field theory 5 is a good description 




(a) Bond centered (b) Site centered 




(c) Fibonacci 



Figure 1 : (Color online) (a) The bond-centered Cayley tree, (b) The 
site-centered Cayley tree. In both cases all sites, other than those 
on the boundary, have coordination 3. (c) The "Fibonacci Cayley 
tree" is constructed hierarchically and has some coordination 2 sites. 
Our choice is motivated by the fact that such a cluster has no local 
sublattice imbalance (refer text for details). All clusters are bipartite 
(the dark(red) and light(green) colors show the two sublattices) and 
have no loops. 



of the many body ground state and can reproduce many of its 
features quantitatively. 

Previous studies of this model by Otsuka 20 and Friedman 21 
focussed primarily on ground state properties and excited 
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states were not considered in these studies. More recently Ku- 
mar et. al 22 have significantly extended this analysis to both 
the spin- 1/2 and spin-1 Heisenberg model. We use all these 
studies as useful benchmarks for our own numerical calcula- 
tions. 

The Cayley tree is a lattice without loops, with all sites (ex- 
cept those on the boundary) having a coordination z. The 
number of sites on the boundary is asymptotically the same 
as the number of sites in the bulk. This is a pathology of 
the Cayley tree, in contrast to lattices in Euclidean space. As 
a consequence, different ways of approaching the thermody- 
namic limit may give different results in any problem based 
on the Cayley tree. In particular, when we use three different 
limiting sequences of finite clusters (see sec. II), we find qual- 
itatively different limiting behaviors; a central thrust of the 
paper is to explain why. But it is generally possible in Cayley 
tree problems to define a "bulk" limit if, rather than averaging 
expectations over the entire cluster, one evaluates them only 
within a cluster of some fixed size which is kept at the center 
while the boundaries get farther and farther away. 

In the literature, the "Bethe lattice" is often used to refer 
to the thermodynamic limit of the "Cayley tree", where the 
effect of boundaries is eliminated (see for example Laumann 
et.al 23 which considers random graphs of fixed connectivity). 
For the purpose of this paper, the presence of open boundaries 
will play an important role in the low energy physics and so 
we will simply use the term "Cayley tree" from now on. 

From a theorist's perspective, the Cayley tree achieves 
many simplifications which makes exact solutions possible. 
For example, the Bose Hubbard model on this lattice has been 
recently solved by Semerijan, Tarzia and Zamponi 24 . It has 
also been used as a tool to study the Brinkman-Rice transi- 
tion for the Hubbard model 25 . It has also found applications 
in the treatment of the quantum impurity problem which is at 
the heart of Dynamical Mean Field Theory (DMFT) 26 . There 
has also been interest in the study of dendrimers 27-29 , but it 
does not appear that a spin model has been realized on such a 
topology experimentally. 

The complete absence of loops makes this lattice conducive 
for the Density Matrix Renormalization Group (DMRG) 
algorithm 30 . With the DMRG method we have an explicit 
(yet compact) representation of ground and excited state many 
body wavefunctions which gives us lots of information to un- 
derstand the low energy properties of these systems. In partic- 
ular, reduced density matrices can be used as tools to under- 
stand properties of these states 31 . 

Various versions of the DMRG have been used 20 22 to study 
the spin 1/2 Heisenberg model on the Cayley tree. In addition, 
Lepetit et. al 32 , have considered the fermionic Hubbard model 
at half filling. More recently, Murg et. al 33 have also used 
the Cayley tree to embed quantum chemical systems to study 
them with tree tensor network algorithms (closely related to 
the DMRG). 

The remainder of the paper is divided as follows. In sec- 
tion II we introduce the model and lattices being investigated 
and define a measure of sublattice imbalance associated with 
them. In section HI, we give a brief overview of our imple- 
mentation of the DMRG algorithm applied to generic trees. 



In section IV, we discuss the general properties of the ground 
state and excited states of the bond-centered Cayley tree. In 
section V, we use a variational ansatz given by the single 
mode approximation, in conjunction with an argument from 
first order perturbation theory, to explain the finite size scal- 
ing of the spin gap. Finally, in section VI, we corroborate our 
observations of the ground state properties, by the use of the 
Schwinger Boson mean field theory (SBMFT). 

II. THE MODEL 

We consider the nearest neighbor antiferromagnetic spin 
1/2 Heisenberg model with uniform coupling J, 

In this paper, we use spin rotational symmetry of the Hamil- 
tonian (2.1), to label many body states by \S,S Z ), where S 
refers to the spin of the state and S z is its z component. 

On a bipartite lattice (with sublattices A and B), like the 
Cayley tree, with sites in sublattice A and n b sites in sub- 
lattice B, it is rigorously known 34 that the ground state of the 
Heisenberg Hamiltonian has a net spin S = | nA ~ nB | . 

The first kind of tree we consider is the "bond-centered" 
Cayley tree of the form depicted in Fig. 1(a). The number of 
sites N s for such a cluster is related to the "generation" g by, 

N s (g) = 2 g+1 -2 (2.2) 

Since the bond centered clusters have no "global imbalance" 
i.e. ua — riB, the ground state is a singlet. We have observed 
that the first excited state is a triplet. 

As mentioned before, the notion of "local sublattice imbal- 
ance" will be crucial in understanding the low energy physics. 
For the bond centered cluster, we define a measure of imbal- 
ance (which we refer to as 1^ from here on) by dividing the 
cluster at the central bond into two equal parts. We count the 
excess of one sublattice over the other in one half of the clus- 
ter and multiply by 1/2 for spin 1/2. It can be easily shown 
that Ib(g) is related to the generation g as, 



where +(— ) is for g odd(even). 

Fig. 1(b) is the more usual way of defining a Cayley tree 
and which we refer to as "site-centered". The number of sites 
is related to the generation g by, 

JV a (ff) = 3(2»-l) + l (2.4) 

Unlike the bond centered cluster, a global sublattice imbal- 
ance exists here which leads to a ground state spin of So = 
2 9_1 . We find it convenient to measure the imbalance I s (g) in 
either of the three symmetrical arms of the site centered Cay- 
ley tree. As before, we measure the excess sites of one sublat- 
tice over the other (in one arm) and multiply by 1/2. This def- 
inition is particularly convenient as it gives us I s (g) = h{ff) 
for all g. 
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A recent publication on the study of the Heisenberg model 
on the Cayley tree by Kumar, Ramasesha and Soos 22 consid- 
ers the site-centered clusters. We confirm their results for 
the site-centered case, but interestingly find that the bond- 
centered cluster has significantly different ground and excited 
state properties. We will provide some brief comparisons in 
section IV to illustrate this point. 

How is the situation different if there is no imbalance lo- 
cally? To address this, we introduce the "Fibonacci Cayley 
tree", an example of which is shown in Fig. 1(c). The recipe 
for generating this lattice is to make generation g + 1 by at- 
taching one branch each from generation g and g — 1. If we 
label the number of odd/even sublattice sites by A g and B g , 
then (counting the root as even), we get, 

A g+1 = 1 + Bg + Bg-! (2.5a) 
B g+1 = Ag+Ag^ (2.5b) 

The total number of sites N s at generation g + 1 is, 

N s (g +l) = A g + B g = l + N s (g) + N s {g - 1) (2.6) 

Observe that N s (g) satisfies the Fibonacci recursion (hence 
the name). The size of this lattice grows as t 9 where r is 
the golden ratio r = (1 + V5)/2 ~ 1.618. Also, every 3rd 
generation is unbalanced by 1 (and every other generation is 
both globally and locally balanced). 

In order to have a balanced cluster at every generation, we 
combine two identical generation g Fibonacci constructions 
(as in equation (2.6)), by introducing a bond connecting their 
roots. 



III. DENSITY MATRIX RENORMALIZATION GROUP ON 
GENERIC TREE GRAPHS 

The Density Matrix Renormalization Group (DMRG) is a 
powerful numerical technique for studying many-body sys- 
tems. It was developed by White 30 for one dimensional 
systems to remedy the problems associated with the Nu- 
merical Renormalization Group (NRG) 35 . Otsuka 20 and B. 
Friedman 21 realized that the DMRG is also applicable to the 
Cayley tree. Their procedure was for regular trees and utilized 
an infinite-system DMRG. A recent publication by Kumar et. 
al 22 (for the site-centered lattice) improves upon the scaling of 
previous algorithms, by considering a efficient way to hierar- 
chically construct the lattice. 

Our implementation differs from all of the above, as it al- 
lows us to study the properties of any finite tree (and not nec- 
essarily regular ones, eg. percolation clusters 36 ). We outline 
the details in the remainder of this section. 

We spell out the notation we have used in this section, d is 
the number of degrees of freedom per site. For example, d = 2 
for a spin 1/2 Hamiltonian. z is the coordination number of 
a site (Although our discussion talks in terms of uniform z, 
a generalization to site dependent z is straightforward). M 
will be used to denote the number of retained states on a local 
cluster of sites (or "block"). 



A. Initialization 

Our algorithm starts with the generation of a suitable guess 
wavefunction. As Fig. 2 shows, this is done by performing 
energy based truncations of the Hilbert space on a successive 
hierarchy of clusters of sites (or "blocks"), assuming they are 
completely disconnected from the rest of the system. By "en- 
ergy based truncation" we mean that on each block we retain 
only the M lowest energy states of the local block Hamilto- 
nian, where M is a parameter that determines the accuracy of 
the calculation. This procedure is carried out in parallel, be- 
ginning from the boundary sites of the cluster and terminating 
when one reaches the geometrically central point (called the 
"focal" point or "root"). 

While the end result of our calculations are expected to be 
independent of this initialization, a good choice can greatly 
accelerate convergence. In addition, when targeting an excited 
state (say in a high S z sector), we often introduce a "fake" 
magnetic field in the Hamiltonian to favor retention of those 
states on the blocks that describe the high energy wavefunc- 
tion. 

At the end of the initialization calculation, one has a crude 
description of the full Hamiltonian in terms of the root degree 
of freedom surrounded by z blocks. 




Figure 2: Warm-up step of the DMRG involving multiple indepen- 
dent renormalization procedures on the tree utilizing energy based 
truncation. The "army continues to march in" from all sides till one 
reaches the geometrically central point (often called the "focal point" 
or "root"). Here we show all stages for a given tree. The red dots rep- 
resent renormalized blocks. 



B. Density matrix based truncation 

We will now consider how every "iteration" in the DMRG 
is carried out to systematically approach the ground (or ex- 
cited) state(s) of the system. For this purpose, we require a 
description of the full Hamiltonian in terms of a "site" de- 
gree of freedom (here (t, 4)) and the basis spanned by the 
M z states retained on the z blocks surrounding it. 

At every iteration we use the Lanczos algorithm 37 to solve 
for the lowest energy eigenvector. Treating one of the blocks 
as the "environment" and the remaining z — 1 blocks and the 
"site" collectively as the "system", we obtain the reduced den- 
sity matrix (RDM) of the system from the ground state of the 
superblock by computing; 

PRDM = Tr cnv (IV'GSXV'GSI) (3-1) 
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Each block takes its turn being the "environment" while the 
other blocks together with the "site" act as the system. (The 
order of choosing environments is not very crucial to the final 
result.) 

In addition to the ground state density matrix, we have also 
targeted higher excited states, by performing a state averaging 
of the reduced density matrix, 

avn 

Prdm — ^■ z -> 

where is the positive weight given to the density matrix 
formed from state \ipi). In most cases, we simply used Wi = 1 
for all states we targeted. An advantage of state averaging is 
that it often helps us to prevent (or get us out of) solutions 
which are local (but not global) minima. The reduced den- 
sity matrix is diagonalized and only M states with the highest 
eigenvalues (out of the total dM z states) are kept on the block. 



C. Sweep algorithm 

The root of the tree has no special significance in the algo- 
rithm once the initial guess wavefunction has been created and 
a representation of the block states has been generated. For 
example, once the density matrix based truncations with the 
root as "site" are completed, the algorithm proceeds to con- 
sider a new division of the lattice into "systems" and "envi- 
ronment" by considering a site on the shell one step out from 
the root. Each of the z sites connected to the original root gets 
its turn being the "new root". 




Figure 3: Division of the Cayley tree locally into site, system(s) and 
environment as required by the DMRG. One renormalization step 
consists of combining the site and system(s) into a new system, re- 
taining the states governed by tracing out the environment degrees of 
freedom. 

Once we reach the exterior of the cluster, we "sweep in". 
However, this time (and for all future sweeps) we have an en- 
vironment present, whose states we trace over to guide the 
density matrix based truncation process. This in-out-in sweep 
continues till convergence of the energy is attained. Thus the 
scaling of algorithm is zN s m z ~ 1 d where N s is the number 
of sites in the tree. For a highly symmetrical lattice (such as 
the Cayley tree), one can reduce this computational cost to 
zlog(N s )M z ~ 1 d (however this has not been implemented). 
One complete sweep is defined here as a sweep out (from the 
root to the boundaries) followed by a sweep in (from bound- 
aries to the root) operation. 



D. Computing Expectations 

Just as in the usual one dimensional DMRG, we have an ex- 
plicit representation of the wavefunction in terms of the block 
states and we can efficiently compute various kinds of cor- 
relation functions. For the purpose of this paper, we have 
simply measured the spin-spin correlation functions (s*, • Sj) 
and the matrix element (1| sf |0), where |0) (|1)) refers to 
the ground state singlet (first excited state triplet) and has the 
labels |5 = 0, S z = 0} ( |5 = 1, S z = 1)). Both these func- 
tions are needed for calculating the coefficients occurring in 
the single mode approximation, which will be discussed in 
section V. The latter is computed by targeting the ground and 
excited state in the same DMRG run, so that both states have 
the same block representation. 

We also study the eigenvalues of the reduced density matrix 
(for a particular division of the lattice), collectively known 
as the "entanglement spectrum". These turn out to be a very 
useful probe of the low energy degrees of freedom (as we will 
see in section V C). This needs no extra effort in the DMRG, 
since these eigenvalues are computed anyway as part of the 
truncation process. 



E. Parameters and Benchmarks 

All calculations reported here are for the z — 3 case and 
M < 160. For all systems considered here we get reasonable 
convergence of the energy within 20 sweeps. 

To benchmark our calculations we have also compared our 
results with Exact Diagonalization data where possible. In 
particular, for the bond centered tree we considered the ground 
state energies and correlations in all S z sectors of the 30 site 
cluster and some high S z sectors for the 62 site cluster. 

One can see from Table 1 that the convergence of the en- 
ergy and the spin gap is rapid as a function of the number 
of retained states on a block (M) for the site centered (in the 
S z = Sq sector) and Fibonacci trees (in the S z = sector). 
However, for the bond centered case (Table II), the conver- 
gence with M is much slower. We observed that it was only 
towards the M = 100 mark that the SU(2) symmetry of the 
wavefunctions began to be restored in our calculation. While 
larger M calculations are certainly desirable, we believe the 
present calculations are reasonable and confirm the existence 
of an anomalous energy scale in the many body spectrum of 
the bond centered tree. 



M 


N s = 


190 


N s = 


176 




Egs 


A 


Egs 


A 


20 


-74.54049387 


0.9397293 


-76.46983049 


0.08834668 


60 


-74.54054021 


0.9397198 


-76.47071767 


0.08725531 


100 


-74.54054022 


0.9397198 


-76.47072851 


0.08725207 



Table I: Energy (Egs) and Spin gap (A) for the 190 site site-centered 
and 176 site Fibonacci lattices for various values of M, 
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The question of finite size scaling can be considered tricky 
given the small magnitude of the singlet-triplet gap. However, 
by computing the ground state energies in every S z sector, 
and fitting the "tower of states", we believe the errors are mit- 
igated. 



M 


Egs 


A 


<l|a? |0> 


<1K|0> 


60 


-49.3405119 


3.4 x 10" 


-4 


0.310 


0.370 


80 


-49.3412938 


5.7 x 10" 


-4 


0.300 


0.344 


100 


-49.3415002 


6.0 x 10" 


-4 


0.280 


0.337 


120 


-49.3415347 


6.0 x 10" 


-4 


0.279 


0.335 


140 


-49.3415521 


6.0 x 10" 


-4 


0.278 


0.335 



Table II: Energy {Egs), Spin gap (A) and s + matrix elements for 
the central (c) and boundary (6) sites for the 126 site bond-centered 
for various values of M.|0) and |1) refer to the lowest singlet and 
triplet respectively. 



IV. GROUND AND EXCITED STATES 

Using the DMRG algorithm presented in section III, we 
calculate the ground state energy and spin gap, which serve 
to point out the differences between the bond and site cen- 
tered clusters. To highlight the role of local imbalance, we 
also compare the excited states of the bond-centered and Fi- 
bonacci trees, both of which are globally balanced. 

A. Ground State Energy, Spin-spin correlations and Spin Gap 
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Figure 4: Ground state energy per site for the bond centered and 
site centered Cayley trees for various lattice sizes. The Fibonacci 
Cayley tree energies are out of the scale considered. A fit to the 
bond centered results given by Eq.(4.2) has been shown. Inset: Finite 
size scaling of the energy gap A plotted on a log-log scale. The 
bond-centered and Fibonacci clusters appear gapless in the infinite 
lattice limit, with a finite size scaling exponent of a w 2 and a fn 
0.6. However, the site-centered clusters have a finite A in the infinite 
lattice limit in concordance with the results of Ref. 22. The lines 
shown are fits to the DMRG data using equations(4.3, 4.4). 

Despite the dissimilarities between the three lattices, the 
"bulk limit" of the estimated energy per bond, based on tak- 
ing an average of nearest neighbor (s{ ■ Sj) over the innermost 
bonds, is roughly identical for all three kinds of Cayley tree, 
about -0.35J. 



We consider the bond-centered, site-centered and Fibonacci 
clusters and compute the ground state energy and the spin gap 
A defined as, 



A = E GS (S Q + l)-E GS (So) 



(4.1) 



where Sq is the spin of the ground state. 

In order to compute the energy per site in the infinite lattice 
limit eo we fit to the functional form, 



Egs 



e Q + 



b_ 



(4.2) 



Depending on whether we fit the bond or site-centered clus- 
ters we get e = -0.393855(2) and e = -0.393854(2) re- 
spectively both of which are consistent within error bars of 
extrapolation and with the value of cq — —0.39385 reported 
by Ref. 22. 

In comparison, as Table 111 shows, the energy per site of the 
Fibonacci tree is significantly lower than either of the bond 
or site centered trees. This is done by the formation of very 
strong nearest neighbor dimers (especially on the boundary 
as the degree of dimerization dies down as one proceeds in- 
wards). The degree of boundary dimerization is more limited 
in the site and bond centered trees. 



Cluster 


eo 


a 


b 


A 


Bond centered 


-0.393855(2) 


0.292(1) 


-1.05(4) 




Site centered 


-0.393854(1) 


0.291(1) 


+0.099(3) 


0.73 + 2.86/VK 


Fibonacci 


-0.435433(1) 


0.167(1) 


-0.42(4) 





Table III: Ground state energy per site eo, finite size scaling param- 
eters for the ground state energy a, b (from Eq.(4.2)) and spin gap 
A (from equations(4.3,4.4) for the bond-centered, site-centered and 
Fibonacci clusters. 



We now turn to a discussion of the spin gap. The site-centered 
cluster has a spin gap in the infinite lattice limit, which we 
obtained by fitting to, 



A = A n 



c 



(4.3) 



and found A^ to be 0.731(4) and a ~ 0.5. 

The bond-centered and Fibonacci clusters appear to be gap- 
less in the infinite lattice limit, based on cluster sizes up to 254 
and 464 sites respectively. We computed the finite size scaling 
of this gap using, 



A - n: 



(4.4) 



Empirically, the value of a ~ 0.6 for Fibonacci and a ~ 2 for 
the bond centered clusters matches rather well with our data 
(see inset of Fig. 4). 

Denoting the ground state as |0), we also compute the 
ground state connected correlation function defined here as, 

Gy = (0 \3i ■ I 0) - (0 \sf | 0) (0 | S ) | 0) (4.5) 

For the bond-centered and Fibonacci clusters, (0 \ s* \ 0) = 
for all i, and hence we have, 



G* 



(0| 



o) 



(4.6) 



Fig. 5 shows some sample correlation functions on all three 
lattices. The marked difference in the behavior of the spin gap 
and the spin correlations between the site- and bond-centered 
clusters can be attributed to a different manifestation on the 
two type of clusters of the spontaneous spin symmetry break- 
ing occurring in the thermodynamic limit. First, the behavior 
of the spin correlations can be understood in the following 
way: on the site centered clusters, the system has an extensive 
total spin S = 2 9 ~ 1 in the ground state. By choosing a partic- 
ular state out of this large multiplet it is possible to orient the 
Neel vector at no energy cost in this finite size system. In par- 
ticular if one considers the S z = S state of the multiplet, the 
local (sf) expectation values will reflect the ordering pattern 
directly. This situation is somewhat analogous to ferrimag- 
netic systems. In the case of the bond-centered clusters we 
have a unique 5 = ground state, which forbids finite (sf ) 
expectation values on a finite system, and the long-range or- 
der then has to be coded in the correlation functions leveling 
off to a finite value at long distances. 



B. Low energy Tower of States 

For the balanced Heisenberg antiferromagnet (ua = n.B 
with a singlet ground state) on a regular unfrustrated lattice 
(e.g. square in 2D or cubic in 3D), with number of sites N s , 
it has been noted and argued 17 ' 18 that the low energy Hamilto- 
nian can be described by a rigid rotor model, 
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S(S + 1) 
21 



where S it the total angular momentum (spin), and 



(4.7) 



(4.8) 



where \ is the susceptibility of S to a field coupling to it. 
Thus, we have a sequence of multiplets, popularly referred to 
as the "Anderson Tower of States", which become degenerate 
in the limit N s — > oo thus making SU (2) symmetry breaking 
possible in this limit. 

To see if this behavior is found on the Cayley tree, we com- 
puted the ground state energy in every S z sector. This may 
be identified with the multiplet energy E(S), since E(S) is 
mono tonic in S and S > S z . 

For the bond centered clusters, a tower of states exists, but 
an anomalous one. In Fig. 6 we observe that in the bond- 
centered case, the E(S) curve consists of two linear pieces, 




(b) i 



o 

(3 



0.31 
-0.10 
3.03 
0.01 



4 5 6 7 8 9 10 11 12 
Distance from tip spin 







, . . , 2V 
.... \ 
. *. .* . 

• . x « • 

w 


H G «..'I 



2 3 4 5 

Distance from central spin 



0.4 
0.32 
0.24 _ 
0.16 - 
0.082 
0.002 




4 5 6 7 8 9 
Distance from tip spin 



11 12 



Figure 5: (Color online) Ground state spin-spin correlations G a ,i, 
as in equation (4.5), for the three kinds of Cayley tree. One (refer- 
ence) spin a is held fixed while the other spin i is take at distance a 
along the highlighted path. In (a) and (c), the reference spin is at a tip 
(a — > tip), and DMRG results are compared with numerical solutions 
of Schwinger Boson mean field theory (SBMFT) from Section VI. 
In (b), the reference spin is at the central ("root") site (a — > 0). (a) 
Bond-centered tree with N s = 126 sites. The SBMFT corrrelations 
shown have been scaled up by an overall factor of 1.8 to take into ac- 
count corrections beyond mean field (in the broken symmetry phase). 
The DMRG and SBMFT results show good agreement, asymptoting 
to a constant, (b) Site-centered cluster with N s — 190 sites, in the 
maximum S z member of the ground state multiplet (S z = So). The 
magnetization | (sf) | is also shown, as a function of distance from 
the center. Even though the connected correlation function decays 
to zero exponentially fast, the long range order is encoded in the 
fact that that the magnetization is non zero, (c) Fibonacci tree with 
7V S = 40 sites. For the "quantum disordered" phase, the SBMFT 
correlations were scaled up by an overall factor of 3/2 (For details 
see sec. VI B). Correlations appear to be decaying exponentially with 
distance. 
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Figure 6: Lowest energy level in every S z sector for the 108 Fi- 
bonacci and the 126 site bond-centered Cay ley tree. The range of S 
from to S* has been magnified and shown in the inset for the 126 
site cluster. It shows a tower of states with a much larger moment 
of inertia than expected from the Anderson tower. This is seen as a 
sharp kink like feature at S* . In contrast, for the Fibonacci tree, the 
transition from the low to high S behavior is less well defined. 



joined at a critical spin S* which depends on the cluster size. 
In effect, the system has two moments of inertia, Ji ow for S < 
S* and the (much smaller) / h i g h for S > S*. Fig. 7 plots 
the size dependence of these moments of inertia, showing that 
/low ~ while /high ~ N s ; it will be our task in Sec. V C, 
to explain this difference. We also observe that 

S* = 2I b (4.9) 



where J5 is the sublattice imbalance on one half of the bond- 
centered cluster as defined in equation (2.3). 
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Figure 7: Moment of inertia of the low (Ji ow ) and high energy (/high) 
rotors as a function of lattice size (N s ) for the bond centered tree of 
various sizes as shown on a log-log plot. 

For the Fibonacci tree, we do not find a clear distinction 
between the two linear pieces as seen in Fig. 6. The scaling of 



the energy gaps (E(S + 1) — E(S)) changes from N s 6 for 
small S, to 1/N S Anderson scaling for large S (4.7). 

In contrast with the gapless spectrum of the bond-centered 
and Fibonacci clusters, the site-centered case has a finite spin 
gap (see Table III) in the infinite lattice limit. 

A complementary probe of the low energy physics is the 
magnetization (m) (see Fig. 8) defined as, 

™ = ^X>*> Gs (4 - 10) 

i 

as a function of a uniform applied magnetic field h. We ob- 
serve a rapid increase in magnetization for small h and it 
seems the saturation magnetization is about m* ~ 1/6 (i.e. 
m* /m sat ~ 1/3). Beyond this rapid rise of the magnetiza- 
tion at small field, the magnetization curve displays a sur- 
prisingly rich structure, with several plateau-like features at 
intermediate magnetization, linked by more susceptible parts 
of the magnetization curve. We note that the first plateau at 
m* seems to have a similar extent in magnetic field as for the 
site-centered clusters studied in Ref. 22. The detailed charac- 
terization of the magnetization curve as N s — > 00 appears to 
be an interesting problem for future studies. 



E 




t , 1 , 1 , 1 , 1 , 1 , 1 

0.5 1 1.5 2 2.5 3 

h/J 

Figure 8: Magnetization curves for bond-centered Cay ley trees of 
various sizes obtained using DMRG. Inset: The rapid rise of the 
magnetization with application of a small magnetic field indicates 
a susceptibility diverging with system size. 



V. SINGLE MODE APPROXIMATION FOR THE EXCITED 

STATE 

Can we understand the origin of these "anomalous" states 
within a unified framework for treating sublattice imbalance? 
In order to answer this, we study the nature of the triplet ex- 
cited state and its relation to the ground state using the sin- 
gle mode approximation 38 (SMA for short). Assuming the 
knowledge of the ground state wavefunction |0) (analytically 
or from a numerical method), the SMA ansatz for the trial 
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12 14 16 18 20 22 24 26 28 30 

Total S z 

Figure 9: Magnetization curves for sites on various shells on the 62 
site bond-centered Cayley tree. The subscript refers to the shell on 
which the site is present, that is, refers to the central two sites and 
the 4 refers to the boundary. The boundary spins have a high local 
susceptibility. 



state |1') = \S = 1, S z = 1) state is given by, 



II') = 



— J2u r s+\0) (5.1) 

3=1 



where Uj are variational parameters and Afy is a normaliza- 
tion factor given by, 



Afy 



'Y, u *k ui (°\ 



o) 



(5.2) 



Using the spin symmetry of the Hamiltonian (and hence its 
eigenfunctions) and the fact that the ground state has S z = 0, 
the normalization factor Afy (5.2) can be written as, 



(5.3) 



fc.i 



where is the spin-spin correlation function previously de- 
fined in Eq.(4.6). 

For a singlet ground state, observe that, there is a gauge 
degree of freedom in the choice of the SMA wavefunction 
(5.1), an arbitrary constant shift u, i.e., 

5> iS +|0> =5^ + ^10) (5.4) 



since for a ground state with total spin S = 0, 

«JXlQ>=«fi&|o> = o 



(5.5) 



It can also be shown 51 that the normalization Afy in Eq.(5.3) 
is invariant under the transformation Ui — > + u. 

For a given trial state ?/'t which is a function of some pa- 
rameters Uj, the variational principle guarantees that, 



En 



<V>t({M)I<MK})> 



(5.6) 



where Eq refers to the energy of the lowest lying state with 
the same symmetry as the trial wavefunction. The best wave- 
function is obtained by optimization of the parameters p by 
minimizing the variational energy Et- Note that within the 
SMA formalism, the ground state (and hence the ground state 
energy) is assumed, which implies that the SMA spin gap is 
also a variational estimate for the true spin gap. 

Here we will show that the SMA does turn out to be a 
very good ansatz for this system based on the close to 100 
% overlap of the SMA wavefunction with the wavefunction 
from DMRG. Our procedure does not require explicit knowl- 
edge of the wavefunction, rather only certain matrix elements 
and correlation functions are necessary. Taking intuition from 
numerical calculations, we construct coefficients Uj occurring 
in Eq.(5.1) to obtain a simple variational state with a gap that 
goes to faster than 1/N S . The aim of this section is thus to 
understand the operator that creates the triplet wavefunction 
from the ground state singlet. This in turn will tell us how 
the spins collectively act, accurately predicting the existence 
of an "anomalous" energy scale. 

A similar SMA calculation was performed by Laumann et. 
al 39 , who considered AKLT models on the Cayley tree. In- 
stead of making their assumption of exponentially decaying 
spin-spin correlations, we use the values of (1| sj~ |0) and Gij 
from our DMRG calculations. In addition, we make no as- 
sumptions about the variational parameters {uj}. 



A. Obtaining the SMA coefficients from maximization of 
overlap with the true wavefunction 

The overlap of an approximate wavefunction with the ex- 
act one can serve as a good measure of its quality. Thus we 
consider the overlap of the SMA wavefunction with the true 
triplet state |1) = \S — 1, S z = 1) i.e., 



O 



fl|l M 



Efe,; lu k uiG 



ki 



and try to maximize it. We have defined fa to be, 

/i = <i|4l°) 



(5.7) 



(5.8) 



We christen this term the "flippability" of a spin. This is mo- 
tivated from its very definition: the more easily flipped spins 
have a larger contribution to the formation of the first excited 
triplet. 

Now we present a method to obtain the optimal parameters 
Uj to construct To meet the requirements of a high over- 
lap of the SMA wavefunction |1') with the exact many body 
triplet |1) , subject to the constraint that it is normalized, we 
devise a cost function Csma defined as, 



C s 



MA 



•5>i/* + A(.A/?,-l) 



(5.9) 



where we have introduced A as a Lagrange multiplier. Taking 
{uj} as our variational parameters, and setting the derivatives 
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of Csma to 0, we get, 

dCsMA 



dui 



fi + ^A^ujGa + Xui = (5.10) 



Thus we get the set of equations (one equation for each i) , 



(5.11) 



To explicitly obtain a state |1') which has a high overlap with 
|1) , we solve the above linear equations for u, numerically. 
Note that the matrix G always has exactly one zero eigenvalue 
because of the gauge degree of freedom (5.4). Hence we can 
not simply invert G to obtain Uj. Instead of inverting G, we di- 
rectly solve the linear system of equations (5.11) using dgesv 
inLAPACK. 

A natural choice of gauge for the parameters {uj} is to sat- 
isfy, 



E" 







(5.12) 



We consider various functional forms for Ui and numerically 
compute their overlap with the exact triplet and compare SMA 
gap estimates. 

For the nearest neighbor Heisenberg model, the SMA gap 
is given by (for a derivation refer to Appendix A), 



SMA 



E<fc,i) (uk-uif 
-Y.,, UiUjGij 



G 



kl 



(5.13) 



Observe that the numerator and denominator (being propor- 
tional to Afy) are invariant under the transformation Uj — > 
Ui + u. 

To minimize the SMA gap, one would like to minimize the 
numerator and maximize the denominator of equation (5.13) 
(note both the numerator and denominator are positive). To 
minimize the numerator, we can try to keep sa ui for as 
many bonds as possible, and hence consider the "left-right" 
ansatz, 

J+l i G left of central bond 
1—1 i 6 right of central bond 




Joundary 



Center 



Shell 



Figure 10: (Color online) Amplitude of the SMA coefficients m for 
various shells (normalized with respect to amplitude on the bound- 
ary) for the N 3 = 30,62 and 126 site bond-centered lattices. Inset: 
The sign structure of m is the same as equation (5.14). Dark(red) 
and light(green) indicates negative and positive Ui respectively. 

Our observation from the numerical solution of equation 
(5.11) for the bond-centered Cayley tree is that m > for 
i on one side of the central bond and m < on the other 
side. We have also plotted the amplitudes of the optimal SMA 
coefficients for the 30, 62 and 126 site balanced Cayley trees 
in Fig. 10. 



B. Comparison of various SMA wavefunctions 

We next try to understand the qualitative nature of the SMA 
solution from the perspective of minimizing the triplet energy. 



which is consistent with the gauge condition (5.12). This 
sign structure is in concordance with the numerical solution 
of equations (5. 1 1) for the bond centered Cayley tree. 

Note that this is quite contrary to the "staggered pattern", 
one obtains by solving equations (5.11) for the square lattice. 
The staggered pattern is defined as, 

{+1 i G even sublattice 
— 1 i G odd sublattice 

The staggered pattern is an energetically expensive solution 
for the bond-centered Cayley tree. Even though it maximizes 
the denominator making it O(N^), the numerator is also large 
i.e. 0(N S ). Thus the SMA gap scales only as 0(l/N s ) 

Table IV verifies the arguments presented above by explic- 
itly listing out the SMA gap and overlap with the exact wave- 
function for the various choices of Ui we have considered here. 
Our inference is that the optimal and the "left-right" ansatz are 
qualitatively similar and yield a much smaller SMA gap than 
the "staggered" ansatz. 

The SMA calculations suggest that all the spins are in- 
volved in the construction of the first excited state from the 
ground state. The central spins contribute roughly a third as 
much in the "optimal solution" and exactly as much as the 
boundary spins in the "left-right" ansatz. The point to note 
here is that in either case the contribution of the spins in the 
interior is not small. This suggests that the antiferromagnetic 
bonds between successive shells do a reasonable job of lock- 
ing spins together (within each half of the bond centered tree), 
resulting in an emergent degree of freedom which is what we 
call a "giant spin". This interpretation will be established next 
in section V C. 
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Table IV: SMA gap and wavefunction overlap with excited state from 
DMRG for various functional forms of m 
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Figure 1 1 : (Color online) The entanglement spectrum for the bond- 
centered,site-centered and Fibonacci trees as shown in the inset. A 
refers to the eigenvalue of the reduced density matrix of the shaded 
area. The ground state for the bond-centered and Fibonacci clusters 
was a singlet and only the S z > sectors are shown. For the site- 
centered cluster we chose to work with the maximal S z sector (which 
is the S z =16 for the 94 site cluster), h and I 3 denote the "imbal- 
ance" metric defined in the introduction (refer to Eq.(2.3)). For the 
bond-centered case, the largest degenerate eigenvalues of the reduced 
density matrix indicate a multiplet whose spin length exactly equals 
the imbalance For the site-centered case, the density matrix has 
largest weight in a state whose spin is I 3 . For the Fibonacci case, a 
spin 1/2 state has the largest weight in the density matrix. 



C. The "Giant Spins" Picture 

As we inferred previously, there are indications that strong 
antiferromagnetic correlations force all spins in one half of the 
bond-centered cluster to act collectively as a single magnetic 
moment. We make this understanding more concrete in the 
present section. 

We divide the bond-centered Cayley tree into two equal 
halves at the central bond. Using the ground state, we com- 
pute the reduced density matrix prdm (see Eq.(3.1)) of one 
of the halves and diagonalize it. The corresponding eigenval- 
ues are arranged by total S z and the resultant plot is the "en- 
tanglement spectrum". The calculation of the entanglement 
spectrum needs no extra effort in the DMRG procedure we 
outlined in section III. Appropriate cuts are also chosen for 
the site-centered and Fibonacci trees as shown in Fig. 1 1 . 

The entanglement spectrum shows a copy of the largest 
eigenvalue in every S z sector ranging from —If, to where 
lb is the net sublattice imbalance and is given by (2 9 ± l)/6 
as mentioned in equation (2.3). This indicates the presence of 
a "giant spin" of spin length /{, whose multiplet is given by the 
eigenvectors corresponding to these large eigenvalues. Given 
this picture, we explain the existence of the "anomalous" scale 
of energies. 



T L Tr 




(a) Bond centered (b) Site centered (c) Fibonacci 



Figure 12: (Color online) The "giant spins" are the low energy de- 
grees of freedom for the bond and site centered clusters that arise out 
of sublattice imbalance. 
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1. Bond-centered tree 

The Heisenberg Hamiltonian on the bond centered Cayley 
tree may be rewritten as, 

H = Hleft + H r ig ht + Js ■ Si (5.14) 

We treat the term corresponding to the central bond Js~o ■ si as 
a perturbation within degenerate first order perturbation the- 
ory. The many body ground state on each half is a degenerate 
multiplet of spin I b . Since all spins on the left and right con- 
tribute in a definite proportion to the "giant spin" operators 
Tr, and Tr, one can re-express the expectation values of the 
J so ■ s*i in terms of Tl and Tr. Note that expectation values 
of the term Jsq ■ s\ are computed in the product basis of the 
two systems given by \T L , T[ ) ® \T R , T R ) 

If all the spins in one half of the cluster had equal participa- 
tion in their collective multiplet (observed in the entanglement 
spectrum),then, 

(t l , T[ \s a \ T L T' L Z ) = - (t l , Ti \s e \T L T L z ) (5.15) 

where e (o) refers to any even(odd) sublattice site in the same 
half of the bond centered cluster. Therefore, if one were to 
create the equally weighted spin operator Tl = J2i *i anc ^ 
consider its matrix elements one would get, 

(Ti) = ±2I b (s ) (5.16) 

where refers to the central site in one half of the bond cen- 
tered cluster. The sign depends on whether the central site 
and the boundary sites are on the same (+) or opposite (— ) 
sublattices. 

What happens when the spins are not equally participating 
in the multiplet? The simplifying assumption we make here is 
that within the projected low energy subspace, each individual 
spin half operator at lattice site i, Si, is proportional to the 
operator Tt/ij,. Using the fact that If, ~ N s , this relation can 
be expressed as, 

So = f T L (5.17) 

where the constant 7b has been used to denote the proportion- 
ality factor. A similar relation exists for si and Tr. From 
equation (5.17) it follows that, 

J{s -s 1 )= J ^{T L -T R ) (5.18) 

The Hamiltonian T^ • Tr is diagonalized by coupling the 
left and right "giant spins" into a spin for the whole system 
i.e. T = Ti + Tr, whose eigenstates are given by \T, T z ). 

(T l .Tr) |t ^ = l - (T 2 - T 2 L - T 2 R ) (5.19) 

= ~T(T + l)-I b (I b + l) (5.20) 

where T varies from 0, 1 , 2I b . Note that Tj, and Tr are 

constant and equal to I b . The term Ib(h + 1) is a harmless 



constant energy shift to all states. Thus, the energy spectrum 
(up to a overall shift) as a function of T is simply, 

^( T ) = ^T(T+1) (5.21) 

This is exactly the Hamiltonian of a quantum rotor with a 
"anomalous" moment of inertia scaling as N^. This simple 
picture, hence, rather remarkably explains our numerical ob- 
servations in section IV of the paper. We find tj to be ~ 3.24 
from fits to our numerical data based on Fig. 7. 

Note that though the giant spins are interacting via a weak 
bond, the fact that they are paired up in a singlet ground state 
makes the state highly entangled. In comparison, as we shall 
soon see, the ground state of the site-centered tree (in the 
S z = So sector) is closer to a product state of the giant spins 
and hence has a lower degree of entanglement. (This also ex- 
plains why the convergence of DMRG calculations is more 
rapid with the number of states M for the site-centered case 
as compared to the bond-centered case.) 



2. Site-centered tree 

Let us now perform essentially the same analysis for the 
site centered Cayley tree to further shed light on (and vali- 
date!) the "giant spins" picture. Rewriting the Heisenberg 
Hamiltonian as, 

H = H 1 +H 2 +H 3 + Js • (Si + s 2 + s 3 ) (5.22) 

As before, make the substitution of si,S2,S3 in terms of the 
giant spins Ti,Ta,T3 each of which has spin length I s . Then 
couple the three giant spins into a bigger giant spin T. The 
Hamiltonian now reads as, 

H = H 1 + H 2 + H 3 + ^s ■ T (5.23) 

where 7 S is a constant and T = Ti + T 2 + T 3 . Note that the 
angular momentum coupling rules predict that the value of T 
are in the range from 0,1... to 3/ s . 

Let us now couple the giant spin T to the spin 1/2 degree 
of freedom at the center of the cluster. The energy (in units of 
J"f s /N s and up to a constant of — f ) for the S = 31 s + \ and 
3/ s — | states is given by, 




- 3J S (31. + 1) 

(5.24a) 

- 3J S (31, + 1) 



(5.24b) 

where the superscript indicates the global value of the spin S 
and the subscript is used to indicate the T value it was made 
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of. Similarly we have, 



E. 



31,-1 



37 s -~J (••!/- • ) -3i,(.3/„ 1) 



(5.25a) 



2 

3/ a -l 



(5.25b) 

Simplifying equations and incorporating the constant of — § 
we get, 



3/a + | 

3J S 

31. -J 
3/ s 

3/ s -l 

,3I.-| 
'3I.-1 



3J„ 



-3/ s 



31, 



-31, 



(5.26a) 
(5.26b) 
(5.26c) 
(5.26d) 



Observe that E. 



,3/,-i 
3/ 8 



is the lowest energy. This is in con- 



cordance with the Lieb-Mattis theorem 34 i.e. the ground state 
has total spin So = 3I S — \ . The energy gap (now in absolute 
units) of the So to So + 1 transition is given by, 



A (S -> S + 1) = ^ (61, + 1) 



(5.27) 



of, 



Since I s is approximately N s /18 for large N s we get a gap 



A (S -> S + 1) 



J Is 

3 



(5.28) 



This is consistent with our numerical observation that the 
gap is finite in the large N s limit. Since the measured gap is 
~ 0.73 J we infer that 7 must be ~ 2.19. 

We give further credibility to our giant spin interpretation 
by testing the prediction of the gap for the So to So — 1 tran- 
sition. This turns out to be gapless in the large N s limit, 



A (S 



1) 



J Is 



(5.29) 



which is consistent with our DMRG calculations as shown in 
Fig. 13 . The measured 7^ from the fit of the DMRG data 
to A = Jj s /N s is found to be ~ 2.19 consistent with the 
estimate from equation (5.28), serving as another check of the 
theory. 



3. Fibonacci Cayley tree 

The entanglement spectrum of the Fibonacci Cayley tree 
(see Fig. 11(c)) indicates the creation of a spin 1/2 degree of 
freedom as opposed to the "giant spins" encountered previ- 
ously. The cut shown in Fig. 1 1(c) shows a region having an 
imbalance of one, which is the maximum possible for any cut. 

We believe the lowest energy excitation of this system in- 
volves a breaking of a dimer and creation of two unpaired 
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Figure 13: (Color online) For the site centered clusters, the scaling 
of the So to So — 1 energy gap serves as one of the test beds for the 
"giant spins" picture presented in the text. Here we see a reasonable 
agreement of the DMRG data with the predicted scaling of 1/N S 



spins. Since the bonds in the interior have a decreasing dimer- 
ization strength, the energy required to create this excitation 
is expected to be vanishingly small in the infinite lattice limit. 
This is seen in the spin gap in Table III, but an explanation of 
the observed numerical exponent is a subject of further inves- 
tigation. 



VI. SCHWINGER BOSON MEAN FIELD THEORY FOR 
SINGLET GROUND STATES 

Can we understand the presence or absence of long range 
order on these trees at the mean field level? For this we appeal 
to the Schwinger Boson Mean Field Theory 5 which is capable 
of describing quantum disordered and ordered states within 
the same framework 440 . In this section we will see that this 
theory is a good description of the singlet ground states of the 
bond centered and Fibonacci trees. This section also serves 
to expand the domain of application of the Schwinger Boson 
formalism to situations where multiple parameters need to be 
optimized simultaneously 41 (such as non uniform systems or 
systems with very few symmetries). 



A. Notation and formal set-up 

The SU (N) Heisenberg Hamiltonian is expressed in terms 
of Schwinger Bosons by defining a bond operator, 



N 



(6.1) 



where each Schwinger boson operator 6j m carries two labels, 
i or j for site indices and m for flavor index. The physical 
Heisenberg model Eq.(2.1) corresponds to N = 2. The pro- 
cedure is to decouple the quartic bosonic Hamiltonian into 
a one-body mean field Hamiltonian by doing an expansion 
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in 1/N, Solving the mean field Hamiltonian and putting in 
N = 2 allows us to compare the SBMFT results with DMRG 
calculations. 

The SU (N) Hamiltonian in terms of Schwinger bosons is, 



H 



Heis. 



(a) 



2S' 



(6.2) 



The mapping of the spin Hamiltonian to Schwinger bosons is 
exact if we meet the condition site by site, 



N 
m—1 



NS 



(6.3) 



which ensures that the Hilbert space of the bosons is restricted 
by the spin size S (and the corresponding \S, S z ) states). 
However, we will impose this constraint only on the expec- 



tation b\ m b im 



~^ (^imPim) mf- As a result of not satisfying 
Eq.(6.3) exactly, the mean field energy Em f an d correlations 
differ from exact results (DMRG calculations) by overall scale 
factors 5 ' 42 . 

The decoupled mean field Hamiltonian Hmf is expressed 
in terms of the following variational parameters: a set of 
bond variables Qij for every bond Lagrange multipli- 
ers Xi which impose Eq.(6.3) and a condensate field = 
8i.m(bim) I VN which develops a non-zero value in a phase 
with LRO 40 . Umf is then given by, 6 . 



N 



H 



MF 



E A M Y, b tmbim-NS 
i=X \m— 1 

i<j 

j^^(4>*ihi + 4>ibl 



N 



(6.4) 



The field 4>i couples linearly to Schwinger bosons and is 
conjugate to <5i, m (&] ), As a result of this parametrization, 
the lowest spinon frequency mode ujq of the m — 1 flavor 
develops a macroscopic occupation of Schwinger bosons on 
Bose condensation. At the mean field level the different boson 
flavors decouple and the part of the mean field Hamiltonian 
quadratic inbosonic operators 6j m , b\ m (m = 2, N) can be 
expressed as N — 1 copies of a quadratic Hamiltonian W MF . 
Uf [F is given by, 



u 



MF 



i=i 



imPim 



i<j 



(q 



■6 f 6 f 

1 J im jm 



h.c 



(6.5) 



Integrating out the bosonic fields then gives us a set of 
single-spinon frequency modes and the total mean field en- 
ergy Emf- Since we do not have the luxury of momen- 
tum space, we adopt a real space Bogoliubov diagonalization 
procedure 43 . Since H MF is block diagonal in the flavor basis 



(the Hamiltonian is the same for all N — 1 flavors) we can now 
drop the flavor index m and express it as 



H 



MF 



bt b 



A Q 
Q A 



b 

bt 



(6.6) 



where m ^ 1, b = b%, b^ s ) and A and Q are N s x N s 
matrices given by Ay = XiSij and = Qij for nearest 
neighbor sites i, j. H MF can now be diagonalized by intro- 
ducing Bogoliubov transfomation matrices, U and V defined 
as follows, 



b 

bt 



U V 

v* u* 



a 



a 



(6.7) 



where a = (ai,ct2, ...,aN,) is a vector of Bogoliubov quasi- 
particle annihilation operators. Each quasiparticle creation 
(annihilation) operator at (a p ) creates (destroys) a bosonic 
quasiparticle in the single particle mode /i, where /.i goes from 
1 to N s . The transformation (6.7) allows us to switch to the 
Bogoliubov quasiparticle basis where H MF becomes, 



N ' / 1\ 

u=l ^ ' 



(6.8) 



We can now do a Legendre transformation to replace the field 
cj>i by Pi, allowing us to express Emf as, 



E 



N 



MF 



N 



E 



^ 1 



1 



-LOn . H 

2 J 



i 



(6.9) 



and E cond is given by, 



E cond = e ^iai 2 + E {QnPlPj + h - c -) ( 6 - 10 > 



1=1 



We now consider the case of = 2. The variational pa- 
rameters {Xi}, {Qij} and {&} are determined by minimizing 
Emf with respect to each of them giving the following con- 
straints, 



dE 



MF 



dX % 



^ \/3 t \ 2 + (b\ 2 b l2 ) = S (6.11a) 



dE 



7Y 



dQij 
OEmf 



-Q tl =0 (6.11b) 



i n.n.i 



One of the obvious considerations of applying this theory to 
such a non-uniform lattice is the large number of variational 
parameters which, in general, scale with the system size N s . 
However, due to the symmetries of the Cayley tree the num- 
ber of independent parameters are reduced to order g. The 
task then is to find an optimal set of parameters {A* , Q* 3 , f3* } 
which satisfy the constraints in (6.11a), (6.11b) and (6.11c). 
This is done numerically and is discussed in subsection C. 
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The resulting optimal E M f(K> Q% ■> Pt) can t> e related to the 
physical Heisenberg energy via, 



E 



Heis. 



= 2E 



MF 



JS 2 



(6.12) 



A note on the Pi minimization constraint (6.11c): a trivial 
solution of this equation is to choose pi = for all sites. 
This is the quantum disordered phase. A non zero value of 
Pi signals condensation of Schwinger bosons and long range 
order. 

For finite uniform systems there is no spontaneous sym- 
metry breaking and correspondingly no condensation of 
bosons 44 . Introduction of the condensate field Pi is analo- 
gous to applying a staggered field in the system that couples 
to the staggered magnetization. This breaks the degeneracy 
of the single particle energies corresponding to the two boson 
flavors (for N — 2). The condensate fraction begins to build 
up in the flavor mode with the lowest frequency. 

The algorithm tries to initially find a self-consistent mean 
field solution by varying only the set of A; and Qy's. How- 
ever, if we cannot satisfy the constraints in equations (6.11a, 
6.11b) (with 0i = 0), we resort to adding the condensate field 
Pi, as an additional set of variational parameters. While we 
can not completely rule out the possibility of a solution with 
Pi = 0, we believe that the appearance of a condensate is 
physical. 



B. Correlation functions in Quantum disordered and LRO 
phase 

Here we compute correlation functions that enter into the 
self consistency equations (6.11a, 6.11b, 6.11c). The boson 
density of a given flavor at site i is given by (suppressing the 
flavor index), 



simplicity of notation we define the following combinations 
of U and V matrices, 



Q 

u 

V 



VV T 



(6.15) 



Spin correlations in the Quantum disordered and the broken 
symmetry LRO phase are then given by, 



(Gij) 
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-{6.16b) 



for (i, j) G A or B 



In the quantum disordered phase SBMFT overestimates Gij 
by an overall scale factor 5 of 3/2. We normalize the SBMFT 
correlation function by this factor to take into account the 1 /N 
fluctuation effects for N — 2. Similarly, in the phase with 
LRO we find that we need to scale up Gij by a factor of 1.8 
to make it agree quantitatively with DMRG results. Similar 
overall scale factors have been reported previously 42 . These 
overall scale factors can be suppressed by using Gutzwiller 
projected mean field wavefunctions 45 , which is feasible only 
for small system sizes. Such projected wavefunctions have 
also been shown to give energies in agreement with exact 
results 42 ' 45 . 



(b\bi) 



N a 

p,q=l 
p,q=l 



(V*V J 



(6.13) 



The indices p, q run over all single particle modes and we 
made use of (ct p aJ q ) = S p q which follows since the SBMFT 
many body ground state is the vacuum of Bogoliubov quasi- 
particles. Similarly, 



(bibj) = (UV 1 



(6.14) 



Spin correlation functions G^ can be computed in a similar 
fashion. The only complication arises in the SU(2) broken 
symmetry phase where, due to loss of spin rotational invari- 
ance, we need to compute the (SfSJ) correlations. This in- 
volves evaluating a quartic expectation which we decouple 
into a series of 2 point functions using Wick's theorem. For 



C. Numerical Implementation 

Using the symmetry of the bond-centered Cayley tree, we 
reduce the number of variational mean field parameters. A 
first simplification results from the fact that all sites within a 
given shell on the lattice are equivalent and are therefore as- 
signed the same A,; and Similarly, all bonds connecting two 
successive shells are equivalent and have the same Q'^s. For 
the Fibonacci cluster there are fewer exact symmetries (only 
reflection about the central bond) compared to the Cayley tree 
and therefore a larger number of variational parameters are 
required. 

We use the Nelder Mead Simplex optimizer from the GSL 
library to minimize the following weighted combined cost 
function which aims to reduce (6.9) subject to the constraints 
(6.11a), (6.11b) and (6.11c). Since each of these constraint 
equations are obtained by minimizing (6.9) with respect to 
the variational parameters Aj, Qij and Pi, enforcing the con- 
straints will minimize Emf- 

C ({Xi}, {Qij}, {ft}) = MoCa + ViCq + (6.17) 
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where {fio, fix, /^} are relative weights of terms in the cost 
function and the costs are given by, 

N s 2 

C * = E ((lAI 3 + ( b i2 b i2)) - 5) (6.18a) 
C Q = E + ^(Aij) + jQa) ( 6 - 18b > 

Cp = E + ]T j (6.18c) 

i=i \ j n.n.i y 

In practice, to minimize the weighted cost function (6.17), 
tolerance values for the Cp and C\ are set at 10 -8 and 10~ 14 
respectively and the Qij's are solved for self-consistently. A 
good initial guess for the 's is a pattern that favors dimer- 
ization as suggested by results from Exact Diagonalization for 
small clusters. A good rule of thumb for the bond centered 
cluster is to begin by dimerizing (assigning a high value of 
Qij) the outermost bond and to create a pattern, moving in- 
wards, of alternating bond strengths. 

For cases requiring a larger number of variational parame- 
ters (like in the case of the Fibonacci cluster) it helps to guide 
the Nelder Mead optimization using a "relaxation" algorithm. 
The algorithm starts with an initial guess for the Qi/s and 
allows the optimizer to find an optimal set of A,'s. If the tol- 
erance level for C\ is not met, the Qij's are allowed to relax: 
i.e. the best set of A,;'s is taken as initial guess in an opti- 
mization where the Qij's are now taken to be the variational 
parameters. Thus, by alternating in the choice of variational 
parameters between Aj's and the Qij's for each optimization 
cycle we converge to the desired tolerance limit. 

The stability of the obtained mean field solution was 
checked by adding small random perturbations to the optimal 
mean field parameters. In the quantum disordered phase, the 
saddle point was checked to correspond to a maximum in Aj's 
and a minimum in the Qij's. 

Every optimization cycle scales as ~ N^t, where r is the 
total time for functional evaluation taken by the optimizer to 
converge to a solution. A typical optimization time for the 
N S =126 cluster is about 20 minutes on a personal workstation 
(2.7 GHz Intel Core i7). 



effects (number of sites on the boundary scales as the sys- 
tem size N s ). As a result, for the finite systems considered 
above, the lowest spinon frequency is always gapped ui 7^ 
in spite of a very low (3 cost (~ 10 -8 ). However, with increas- 
ing system size, uq lowers and spinons begin to condense in 
this lowest frequency mode. Fig. 14 shows a finite size fit 
to the lowest spinon mode. Spinon condensation and a very 
small loq suggest long range order in the thermodynamic limit. 



-1.5F 




Figure 14: (Color online) The lowest spinon frequency appears to go 
to as l/iV s .This signals appearance of long range order within the 
Schwinger Boson Theory. 

Since condensation of spinons signals long range order, 
sites with a higher condensate fraction have a greater partic- 
ipation in establishing Neel order on the cluster. By map- 
ping out the condensate fraction on different radial shells in 
Fig. 15(b) we notice the strong correspondence between sites 
with large condensate densities and those with a high "flippa- 
bility" as inEq.(5.8). 

Our results can be put in perspective with respect to a 
heuristic for computing the number of "dangling spins" pro- 
posed by Wang and Sandvik 13 in the context of percolation 
clusters on the square lattice. We "dimerize" the lattice max- 
imally as shown in Fig. 15(a) and the spins that remain are 
called "dangling". These are the representative spins partici- 
pating as the low energy degrees of freedom. (Note that the 
choice of maximal dimer covering is not unique but the num- 
ber of uncovered spins is the same in all such coverings) 



D. Results 

The landscape of the cost function (6.17) in parameter 
space has many local minima with very similar mean field 
energies (differing by 1% or less). As a result, we get a zoo of 
physically plausible mean field solutions, all of which satisfy 
(6.11a) and have comparable Emf- To choose the optimal 
solution from amongst those, we look at the minimization 
constraint (6.1 lc) and hand pick the one which has the lowest 
Cp. In other words, the chosen solution has the lowest spinon 
frequency oj . 

The mean field energy and correlation functions for the 
bond-centered Cayley trees suffer from significant finite size 



VII. CONCLUSION 

In this paper, we have explored the relationship between 
sublattice imbalance and nature of the low energy spectrum of 
the bond-centered, site-centered and Fibonacci Cayley trees. 

For the bond-centered Cayley tree, we find that the spin gap 
scales with size as 1 /Nf where a was found to be s» 2. We 
discover an entire tower of states (Fig. 6) with a much larger 
moment of inertia (Fig. 7) than the Anderson Tower of States. 
This low energy scale persists up to a spin value of S* = 21b, 
where 1^ refers to a measure of the imbalance (or the number 
of "unpaired spins") on the bond centered tree (as in Eq.(2.3)). 
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1 2 3 4 5 
Shell 

Figure 15: (Color online) Above (a): Heuristic for computing the 
number of "dangling spins" as proposed by Wang and Sandvik. Be- 
low (b): The "flippability" as in Eq.(5.8) computed from DMRG and 
the condensate fraction \ fii\ 2 computed from SBMFT on every shell 
of the bond centered Cayley tree. Both quantities are qualitatively 
consistent with each other and confirm the "dangling" spin heuristic 
shown above. 



To highlight the role of sublattice imbalance, we introduced 
the Fibonacci Cayley tree in sec. II, which does not have any 
locally unpaired spins. We found it lacks the low lying states 
characteristic of the bond-centered tree (see Fig. 6). Instead, 
the spin gap vanishes as w 1/7V° 6 . However, both trees have 
similar susceptibilities (~ N s ) at sufficiently large magnetic 
fields. This is because the strength of the dimerization is rel- 
atively weak at sufficiently high energy scales (comparable to 
J), allowing all spins to lock together, leading to an extensive 
susceptibility. 

For the site centered tree, our results are in good agreement 
with a recent study 22 . We report a finite spin gap of A = 
0.731(4) in the infinite lattice limit and a ground state energy 
of e = -0.393854(1). 

Our results can be explained within a unifying framework 
of indvidual spins coupling together to form collective spin 
degrees of freedom which we refer to as "giant spins". The 
idea for coupling big sublattice spins is well known 15 in the 
context of regular lattices. However, we emphasize that the 
"giant spins" are created by coupling all spins (both even and 
odd sublattice spins) in regions with large local sublattice im- 
balances. For the bond- and site-centered lattices, we find 



that all lattice sites have a (nearly) uniform participation in 
the "giant spins". This picture is verified by two independent 
techniques in section V - the single mode approximation for 
the excited state and an observation of the entanglement spec- 
trum. 

In a broader context, our study aims to understand the 
nature of unpaired spins coupling in a strongly dimerized 
background. Such spins, termed as "dangling spins", have 
been predicted from numerical simulations 13 for systems with 
quenched random dilution. Thus, a natural extension of our 
study is to consider such dangling spins on percolation clus- 
ters where we expect that they couple hierarchically to form 
emergent giant spins at different length scales. This will be 
the subject of a future publication 36 . 

In this paper, we have explored several techniques to de- 
velop our understanding. To begin with, we obtain accurate 
many-body wavefunctions and their expectations using an im- 
plementation of the Density Matrix Renormalization Group 
(DMRG) procedure (described in section III) that works for 
generic trees. In general, this procedure is expected to be well 
suited to hierarchical lattices and those with few or no loops, 
such as the Husimi Cactus 46 ' 47 , a lattice of corner sharing tri- 
angles (dual to the Cayley tree and locally like the Kagome 
Lattice). 

To have an understanding at the mean field level, we 
adapted Schwinger Boson Mean Field Theory (SBMFT) to 
a larger number of variational parameters than considered 
previously 42,45 ' 48 . We were able to study spin correlations of 
the bond-centered and Fibonacci trees (singlet ground states). 
Rather remarkably, the theory is quite accurate quantitatively 
in predicting ground state spin-spin correlation functions (see 
Fig. 5), up to overall multiplicative scale factors (as discussed 
in section VI). The recipe outlined in section VI can be used 
to navigate through the zoo of feasible mean field solutions 
by giving relative weights to the constraint equations (6.18a, 
6.18b, 6.18c). 

We believe that most applications of SBMFT have focused 
on quantum disordered phases 49 ' 50 , but the broken symmetry 
phase has received less attention. The setup can also be gen- 
eralized to handle frustrated spin systems without the need 
to have an ansatz for the mean field fluxes or the decoupling 
parameters 49 . This might lead to novel spin liquid ground 
states with new kinds of broken symmetries 41 . 

After writing this paper, we came across recent results by 
G.Misguich , who has done an extensive numerical study of 
SBMFT formalism for spatially inhomogenous states, mostly 
concentrating on gapped systems. However, his study differs 
from ours in that we include the condensate field as a varia- 
tional parameter. It will be interesting to apply our formalism 
to further investigate his proposed set of ground states. 
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Appendix A: Derivation of the SMA gap equation for the 
Heisenberg Model 

In section V, we introduced the single mode approximation 
(SMA) for the excited state in terms of the ground state, 



|l'}=5> Sl +|0) 



(A.l) 



where |0) is the singlet ground state and |1') is the approx- 
imate triplet excited state, i refers to a site index. In this 
Appendix, we will derive an expression for the SMA energy 
gap Asm A m terms of the ground state correlations and the 
parameters Ui. 

The expression for the gap between the ground and first 
excited state is, 



d'li') 
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where E Q is the ground state energy. Plugging in the expres- 
sion for from equation (A.l) we have, 
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Next, consider the nearest neighbor Heisenberg Hamilto- 
nian, 

H = \ E + 2 ( S ^r + ^4)) ( A -6) 

(k,l) v 1 

where (k, I) refer to nearest neighbor pairs. We have included 
a factor of 1/2 outside to compensate for counting each nearest 
neighbor term twice. 



We now calculate [H, sf] occurring in Eq.(A.5). To do so, 
we calculate [sfe ■ si, sf] as, 

[s k ■ s h sf] = s z k [sf,sf] + sf [si, sf] + 

\«i i s T,4]+lbk,4] 4 

= 5 a s z k sf + S lk sisl 

~hi4 s i ~ S *kS z k sf (A.7) 

The numerator of Eq. (A. 5) involves the term (sj [H, sf~\ ). 
Hence we now consider the action of the sj operator on the 
simplified expression for [sfc • s;, sf~\ in Eq.(A.7). Consider 
only the terms that have j = k or j = I (since k ^ I terms do 
not occur in the Hamiltonian we do not have to worry about 
the possibility j = k = I). In addition, time reversal symme- 
try of the ground state wavefunction (equivalent to simply as- 
serting the S z — !• —S z symmetry of the ground state) ensures 
that if both j ^ k and j ^ I then the three point correlation 
function is exactly 0. This latter point is rather subtle and so 
we expand on this in Appendix B. 

Thus the expression for sj [s k • si,sf~\ (after retaining 
only the j = k and j = I terms) is, 
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Inserting (A.8), in the expression for the SMA gap (A.5) for 
the Heisenberg Hamiltonian and utilizing (sf ) = for all i, 
we obtain, 



i-SMA — 
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(A.9) 



Appendix B: Why is (ip\ Sj sf \ip) =0 for distinct j, k, I ? 

To derive the SMA gap equation in Appendix , we used, 

(i>\sjs+sf\ilj) = (B.l) 



for distinct site indices j, k, I. In this Appendix, we will prove 
this statement for any wavefunction which is invariant under 
time reversal. 

Consider three distinct spins i, j, k. Express the wavefunc- 
tion in the basis spanned by the three spins at sites j, k, I and 
the rest of the spins (collectively termed as "environment" e), 
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Since this wavefunction is an eigenstate of the Heisenberg 
model (with no external magnetic fields), it follows that under 
time reversal (denoted by operator T) we have, 



ip — > zip 



(B.3) 



where z is ±1. 

This implies that the coefficients in the wavefunction satisfy 
the relation, 



-sj— sj,— s; 

J — e 



(B.4) 



The action of the operator s sfsf on \ip) from Eq.(B.2) 



yields, 



*74^} = EE w ^i^ 



(B.5) 
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Now acting Eq.(B.5) with from the left and using the 
orthogonality of the basis we get, 
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where e + (e_) reflects the fact that the environment carries a 
net S z of +(— )| since the wavefunction comprises of 5f ot = 
terms only. Under inversion of all spins in e + we get e_. 
With this in mind, consider the second sum on the right. Using 
the time reversal symmetry of the wavefunction i.e. = 

zw\^ and = zw]^} (as seen from Eq.(B.4)), in Eq.(B.7) 
we get, 
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Appendix C: Sch winger Boson Mean Field Theory Calculations 

As mentioned in section VI, optimization of multiple pa- 
rameters occurring in the Schwinger Boson theory for non 
uniform systems was quite a challenging task. Hence, for the 
interested reader, we report the exact values of the parameters 
obtained from our calculations, so that they may be able to 
reproduce our results. 

The optimal mean field parameters are tabulated in Table V 
for different lattice sizes. In each column (from top to down) 
the parameters label inner to outermost most bonds/sites. The 
Qi/s alternate in strength across successive bonds consistent 
with the location of unpaired spins. Similar alternation in the 
condensate field indicates the variation in the density of 
dangling spins across shells. 

The ground state energy from SBMFT for the 126 site clus- 
ter was found to be sa —0.533 J. This is lower than the DMRG 
estimate —0.39385 J. This can be attributed to the well known 
fact 44 48 about the non variational nature of SBMFT energies. 
This is because of not satisfying the constraints in Eq.(6.3) 
exactly. 
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Table V: Optimal SBMFT parameters for bond centered clusters of 
various sizes 



where we have used z = 1. 
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